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A number of studies showed that deep brain stimulation (DBS) can modulate the activity 
in the epileptic brain and that a decrease of seizures can be achieved in "responding" 
patients. In most of these studies, the choice of stimulation parameters is critical to 
obtain desired clinical effects. In particular, the stimulation frequency is a key parameter 
that is difficult to tune. A reason is that our knowledge about the frequency-dependant 
mechanisms according to which DBS indirectly impacts the dynamics of pathological 
neuronal systems located in the neocortex is still limited. We address this issue using both 
computational modeling and intracerebral EEG (iEEG) data. We developed a macroscopic 
(neural mass) model of the thalamocortical network. In line with already-existing models, 
it includes interconnected neocortical pyramidal cells and interneurons, thalamocortical 
cells and reticular neurons. The novelty was to introduce, in the thalamic compartment, 
the biophysical effects of direct stimulation. Regarding clinical data, we used a quite unique 
data set recorded in a patient (drug-resistant epilepsy) with a focal cortical dysplasia 
(FCD). In this patient, DBS strongly reduced the sustained epileptic activity of the FCD 
for low-frequency (LFS, < 2 Hz) and high-frequency stimulation (HFS, > 70 Hz) while 
intermediate-frequency stimulation (IFS, around 50 Hz) had no effect. Signal processing, 
clustering, and optimization techniques allowed us to identify the necessary conditions 
for reproducing, in the model, the observed frequency-dependent stimulation effects. 
Key elements which explain the suppression of epileptic activity in the FCD include: (a) 
feed-forward inhibition and synaptic short-term depression of thalamocortical connections 
at LFS, and (b) inhibition of the thalamic output at HFS. Conversely, modeling results 
indicate that IFS favors thalamic oscillations and entrains epileptic dynamics. 
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INTRODUCTION 

Deep brain stimulation (DBS) for Parkinson's disease (PD) and 
other movement and psychiatric disorders — including dystonia, 
tremor, and depression — is clinically used today as a conventional 
therapeutic procedure for the alleviation of symptoms (Sillay and 
Starr, 2009). Since the early 90s, neurologists also attempted to 
apply DBS to other neurological disorders, typically to intractable 
epilepsies in order to suppress — or at least dramatically reduce — 
the occurrence of seizures [see recent review in Boon et al. 
(2009)]. These studies followed early scientific evidence showing 
potentially beneficial effects of DBS on epileptic neural dynam- 
ics in animal models (Reimer et al., 1967; Hablitz, 1976) as well 
as in patients (Cooper et al., 1973; Davis et al., 1982; Wright and 



Abbreviations: CMN, Centromedian Nucleus; DBS, Deep Brain Stimulation; 
EPSP, Excitatory Postsynaptic Potentials; FCD, Focal Cortical Dysplasia; FFI, Feed- 
Forward Inhibition; GPi, Globus Pallidus; HFS, High Frequency Stimulation; iEEG, 
Intracerebral EEG (depth electrodes); IFS, Intermediate Frequency Stimulation; 
IPSP, Inhibitory Postsynaptic Potentials; LFP, Local Field Potential; LFPsfcd* Local 
Field Potentials recorded in the FCD; LFS, Low Frequency Stimulation; NS, 
No Stimulation; PMC, Premotor cortex; RtN, Reticular thalamic Nucleus, STD, 
Short-Term Depression; STN, Subthalamic Nucleus. 



Weller, 1983). However, contrary to PD, the optimal "antiepilep- 
tic parameters" of DBS for reducing the frequency of seizures are 
much more variable among patients and the number of non- 
responders to stimulation still perplexes scientists. Moreover, in 
responding patients, the fine tuning of stimulation parameters in 
a patient-specific manner remains indispensable for maximizing 
antiepileptic effects. On that account, many fundamental ques- 
tions are frequently raised: where and when to stimulate, at which 
frequency, at which current intensity, and with which current 
waveform? 

The answers to these questions remain bound to our cur- 
rent, and still limited, understanding of the mechanisms by which 
DBS modulates neuronal dynamics, whether normal or patholog- 
ical. Today, the precise mechanisms of neuronal modulation by 
DBS remain elusive. In addition, these mechanisms are controver- 
sial as observed effects are sometimes opposite (Mclntyre et al., 
2004b). Among the many studies reported over the last decade, 
identified mechanisms regarding HFS include: local depolariza- 
tion blockade by HFS (Beurrier et al., 2001), synaptic depres- 
sion due to neurotransmitter depletion (Shen et al., 2003; Kim 
et al., 2012), synaptic inhibition (Filali et al., 2004), disruption of 
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the thalamocortical network's dysrhythmia (Mclntyre and Hahn, 
2010; Kendall et al., 2011). As far as LFS is concerned, some stud- 
ies described a transient synaptic depression that alters synaptic 
transmission (Jiang et al., 2003; Speechley et al., 2007). Finally, 
IFS is routinely used in the context of presurgical evaluation of 
patients with drug resistant epilepsy to map epileptogenic and 
functional brain areas. It has long been observed that this type of 
stimulation is prone to trigger epileptic afterdischarges (Goddard, 
1967). This brief overview shows that the spectrum of involved 
mechanisms is very large and that distinct stimulation frequencies 
trigger distinct cellular/network processes. More precise insights 
into these processes will come with increased knowledge about 
both biophysical and neurophysiological effects of stimulation 
currents on underlying neuronal systems. 

However, the access to cellular and network mechanisms 
induced by DBS is rather difficult in animal models of epilepsy 
and (almost) impossible in patients especially in large-scale sys- 
tems like the thalamocortical loop. An alternative approach is 
the use of computational models based on physiological data to 
first reproduce and then explain changes in cerebral activity as 
a function of stimulation conditions (stimulation site, intensity, 
and frequency). This is precisely the objective of this study, with a 
special focus on the distinct effects of DBS frequency on cortical 
epileptic dynamics. 

Our investigation combines computational modeling and clin- 
ical data. We explored stimulation effects in a lumped-parameter 
mesoscopic neural mass model of the thalamacortical loop, 
inspired from previously published models (Suffczynski et al, 
2004; Lopes Da Silva, 2006; Roberts and Robinson, 2008; Crunelli 
etal.,2011). 

Although these models are lumped representations of under- 
lying neuronal systems, they offer a number of advantages in the 
context of this study. First, neural mass models include subpop- 
ulations of principal excitatory cells and inhibitory interneurons. 
Second, these models were shown to produce realistic activity as 
observed in LFPs or EEG under normal (Freeman, 1973; Lopes 
Da Silva et al., 1974) or epileptic conditions [review in Lytton 
(2008); Wendling (2008)]. Third, main parameters (mean mem- 
brane potential and firing rate) provide access to the investigation 
of several stimulation-induced (patho)physiological mechanisms. 
For instance, a neural mass model was successfully used in the 
context of direct low-intensity pulse stimulation in the hippocam- 
pus to explain the behavior of evoked responses during the 
transition to seizures (Suffczynski et al., 2008). 

In particular, using this model, we analyzed the neurophysio- 
logical effects induced by direct thalamic stimulation on epileptic 
cortical dynamics at low frequency (LF, < 20 Hz), intermediate 
frequency (IF, 20-70 Hz) and high frequency (HF, 70-130 Hz). 
Model parameters were tuned to reproduce a typical pathological 
oscillatory activity observed in a neocortical lesion (focal corti- 
cal dysplasia, or FCD) in a patient with drug-resistant epilepsy. 
Intracerebral EEG (iEEG) signals observed during thalamic stim- 
ulation (centromedian nucleus) of this patient revealed particu- 
larly pronounced frequency-dependent modulation of the FCD 
pathological activity. Therefore, this data set offered the unique 
opportunity to identify key model parameters for which such 
a frequency-dependent modulation could be reproduced and, 



subsequently to get insights regarding the mechanisms under- 
lying the modulatory effects, in the FCD, of thalamic stimu- 
lation. Results revealed that LFS favors feed-forward inhibition 
and short-term depression at the cortical level and that HFS 
inhibits the thalamic activity, while IFS reinforces reticulotha- 
lamic oscillations thus entraining cortical pathological epileptic 
dynamics. 

MATERIALS AND METHODS 

In this section, we present (1) the neurophysiologically- relevant 
computational model that we developed to study thalamic DBS, 
(2) the real depth-EEG dataset used for model tuning and, (3) 
the signal processing methods used for characterizing real and 
simulated EEG signals. 

MODEL OF THE THALAMOCORTICAL LOOP 

In order to study the effects of thalamic DBS on cortical dynamics, 
we implemented a physiologically-plausible mesoscopic model 
of the thalamocortical loop. This model accounts for the aver- 
age activity of both cortical and thalamic compartments which 
include various types of neuronal populations interacting via 
synaptic transmission. This modeling approach was first pro- 
posed in the early 70s (Wilson and Cowan, 1972) and further 
enriched in order to interpret electrophysiological recordings and 
study brain dynamics, in the olfactory (Freeman, 1973) and the 
thalamocortical (Lopes Da Silva et al., 1974) system, for instance, 
as well as the dynamics of cortical oscillations (Nunez, 1974). This 
approach was then developed by other research groups in the con- 
text of state changes in brain dynamics (Wright et al., 1985), visual 
evoked potentials (Jansen et al., 1993), dynamics of the human 
alpha rhythm (Stam et al., 1999) or pathophysiological mecha- 
nisms of ictal transitions in epilepsy (Wendling et al., 2000, 2002; 
Suffczynski et al., 2001; Robinson et al., 2002; Liley and Bojak, 
2005; Breakspear et al., 2006). Later, neural mass models were 
also used in studies dealing with the connectivity among cortical 
regions and the impact of model parameters on the power spec- 
trum of EEG or MEG signals (Robinson et al., 1997; David and 
Friston, 2003; Zavaglia et al, 2006). 

Model architecture 

The model architecture was inspired from previously published 
models of the thalamocortical loop (Suffczynski et al., 2004; Lopes 
Da Silva, 2006; Roberts and Robinson, 2008; Crunelli et al, 201 1). 
In a global view, the model was built of three interconnected 
compartments: a cortical compartment, a thalamic compartment, 
and a reticular compartment, in accordance with previously pub- 
lished models (Figure 1A) and with anatomical data (Figure IB). 
Each compartment includes one or several subpopulation(s) of 
neurons, either excitatory or inhibitory. Generally speaking, the 
input/output functions of a considered subpopulation are rep- 
resented by two mathematical equations that were respectively 
named "pulse-to-wave" (input) and "wave-to-pulse" (output) by 
Walter Freeman (Freeman, 1992). The former is a linear transfer 
function that converts the presynaptic average density of affer- 
ent action potentials into an average postsynaptic membrane 
potential (PSP), either excitatory (EPSP) or inhibitory (IPSP). 
The output function is a static nonlinear function (sigmoid) that 
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FIGURE 1 | Model of the thalamocortical loop. (A) The model 
architecture comprises three main compartments: cortical, thalamic, and 
reticular. The cortical compartment includes three subpopulations: P 
(pyramidal principal neurons), /(" (soma- and proximal-dendrite targeting 
interneurons mediating GABA^ fast currents), and l£, (dendrite-targeting 
interneurons mediating GABA^ s ; ow currents). The thalamic compartment 
represents a generic thalamic nucleus including a subpopulation of 
excitatory thalamocortical (TC) cells. The reticular nucleus (RtN) 



compartment is made up of two GABAergic neuronal populations 
(/{", GABA A fgst currents and /™, GABA^ s / ow ). Excitatory synaptic 
transmission among the considered subpopulations is glutamatergic 
(GLU). (B) Anatomical connectivity of the CMN, PMC, and RtN. This 
diagram represents the anatomy of a particular thalamocortical loop 
interconnecting the CM thalamic nucleus, the PMC, and the RtN. 
Connectivity patterns were inferred from the literature. It is compatible 
with the thalamocortical model diagram presented in (A). 



provides the average pulse density of action potentials fired by 
neurons depending on the sum of EPSPs and IPSPs at the input. 
This non-linear function accounts for threshold and saturation 
effects that take place at the somas and initial axonal segments of 
considered cells. 

Formally, the input function is represented by a second order 
low-pass filter H(s) = W/{s+l/x w ) 2 (where s is the Laplace 
variable). The impulse response of this filter is given by 

W 

h(t) = — ■ t ■ e- t/x » (1) 

X w 

Parameters W and x w are tuned such that h(t) approxi- 
mates the shape of real excitatory (glutamatergic) or inhibitory 
(GABAergic) postsynaptic potentials (Lopes Da Silva et al., 1976). 
The quantity W.x^, is the static gain of filter h. Lumped param- 
eter x w (expressed in s) is linked to the kinetics of synaptic 
currents. It determines both the rise time (f rlse = x w ) and the 
decay time {td eca y = 3.146t w ) of the second order filter impulse 
response h and it is usually adjusted with respect to the phys- 
iological rise and decay times of actual PSPs (Molaee-Ardekani 
et al., 2010). Given the time constantt w , parameter Wean be used 
to adjust the sensitivity of synapses (the maximal PSP amplitude 
is \V.e~ 1 ). An alternative implementation of the h function was 
introduced in Bojak and Liley (2005) and is described in detail 
in Molaee-Ardekani et al. (2013). It is based on a bi-exponential 
pulse-to-wave function with two time constant parameters. This 
implementation allows for the separate adjustment of the rise 
and decay times of PSPs, and therefore a better approximation of 
actual PSPs in some circumstances. Besides, the output function 



is represented by S(v) = — r( e v ° _,,) , where 2en is the maximum fir- 
ing rate, vn is the postsynaptic potential corresponding to a firing 
rate of en and r is the steepness of the sigmoid. 

The cortical compartment 

The cortical compartment was inspired from an existing model 
of the neocortex which proved its capability of generating both 
normal and epileptiform activity. Readers may refer to Molaee- 
Ardekani et al. (2010) for details. In brief, the cortical com- 
partment integrates a subpopulation of pyramidal cells (P, W = 
Ac, x w = x ac in Equation 1) and two inhibitory neuronal popu- 
lations (If and Ij, Figure 1A) representing soma- and proximal- 
dendrite targeting interneurons (GABAa, fast currents, W = 
Gc, x w = x gc in Equation 1) and dendrite-targeting interneu- 
rons (GABA^ s i ow currents, W = Be, x w = x\, c in Equation 1), 
respectively. Pyramidal collateral excitation was implemented as 
in Jansen et al. (1993). 

In addition, these three cortical subpopulations receive exci- 
tatory input from the thalamic compartment. Therefore, feed- 
forward inhibition (FFI) is represented in the model as the two 
subpopulations of interneurons project to the pyramidal subpop- 
ulation (see The Thalamic and Reticular Compartments paragraph 
below). 

Short-term synaptic depression (STD) 

STD is present in the neocortex (Boudreau and Ferster, 2005). It 
can be potentially involved in the context of direct stimulation 
of the thalamus as TC cells directly project to cortical pyrami- 
dal cells. Consequently, this mechanism was implemented at the 
interface of thalamic/cortical compartments. To our knowledge, 
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an implementation of STD in neural mass models has not been 
proposed before. 

In our model, we represented a modulatory effect of the 
amplitude of the average EPSP (parameter Ac) at the level of 
subpopulation P depending on the density of action potentials 
[dAp(t)] coming from the thalamic compartment. This modu- 
latory effect was obtained by multiplying Ac' by a time-varying 
coefficient K(f) e [0.6, 1] where the function describing the evo- 
lution of K(f) was derived from Chung et al. (2002). This study 
shows that: (i) cortical EPSPs drop by 40% under periodic low- 
frequency intense thalamocortical (TC) cell firing and, (ii) this 
drop in cortical EPSP is directly linked to transient depression of 
thalamocortical monosynaptic projections to pyramidal neurons. 

In line with these observations, STD was implemented as fol- 
lows. First, signal d^ p is low-pass filtered (cutoff frequency = 
10 Hz) to restrict the STD effect to LFS. Then, from each time t^ at 

which the filtered signal cf Ap (t) exceeds a firing rate equal to r\, the 
\c(t) coefficient undertakes an exponential decay given by K(f) = 
• e~ f / T where k,, = K(f~) and where t~ is the time instant that 
just precedes t^. The decrease of K.(t) is limited to the time inter- 
val [t^ + 0.45 s] and cannot exceed 40%, total. Parameters r) and 
t were set to 0.8 and 8 s, respectively. 

The thalamic and reticular compartments 

The thalamic compartment was limited to one population of exci- 
tatory neurons (known as glutamatergic thalamocortical - TC 
- cells) receiving glutamatergic EPSPs (W = Aph, t w = x a ph in 
Equation 1) from cortical pyramidal cells (P) and GABAergic 
IPSPs with slow (W = Bjh, x w = Xbrh in Equation 1) and fast 
(W = Gjh> t» = tgT/i m Equation 1) kinetics from the reticular 
compartment (RtN). Here, we increased the time constant (x^jh) 
with respect to x\, c to account for both GABA^ s \ ovi - and GABA^- 
receptor mediated currents in a single variable. TC cells directly 
target both cortical pyramidal cells and interneurons. The acti- 
vation of these GABAergic interneurons subsequently promotes 
inhibition of pyramidal cells after a di- synaptic delay. Therefore, 
TC cells activation induces first an EPSP followed later on by an 
IPSP on cortical pyramidal cells, resulting in feed-forward inhibi- 
tion (FFI). The RtN compartment comprised two inhibitory sub- 
populations, namely if 1 and 7^ T which both receive excitatory 
input from the cortical (W = Ap t , x w = x a p t in Equation 1) and 
the thalamic ( W = Ap t , x w = x a p t in Equation 1) compartments. 

Simulation of stimulation effects 

Stimulation currents induce a perturbation of the membrane 
potential of neurons. At cellular level, this effect can be accounted 
for by the "X£ model", which is well grounded in the bio- 
physics of compartment models (Rattay, 1998; Mclntyre et al., 
2004a; Manola et al, 2005, 2007) (see Miranda et al., 2009 for 
a review) and supported by in vitro experiments (Bikson et al., 
2004; Frohlich and McCormick, 2010). This model AV » X.E 
approximates the membrane potential variation AV as a linear 
function of the electrical field E induced by stimulation (X repre- 
senting the membrane space constant). In our neural mass model, 
the situation is less straightforward as space is not explicitly rep- 
resented, conversely to detailed or mean-field models. However, 
within a certain range of intensity values, it has been shown that 



the membrane potential variation A V is modified in a linear way 
with respect to the electrical field which is itself proportional 
to the stimulation intensity (Bikson et al., 2004). These consid- 
erations led us to also assume a linear variation for the mean 
membrane potential as a function of stimulation intensity, in 
stimulated sub-populations of neurons. In addition, stimulation 
was represented by a train of periodic monophasic depolarizing 
pulses. The pulse width was fixed to 1 ms (as in clinics). Pulses 
were low-pass filtered to account for the average time of repolar- 
ization (set to 4.8 ms) in stimulated sub-populations of cells. The 
resulting stimulation signal was added to the mean membrane 
potential of neuronal sub-populations included in the thalamic 
(TC) and reticular (lf T and I^ T ) compartments of the proposed 
model. The depolarizing effect was weighted by three coefficients 
StCi Sat l and Spa (Table 1 ) accounting for the possibly different 
stimulation impact at the thalamic and reticular level. 

Model parameters, outputs, and implementation 

Parameter values as well as physiological interpretation are pro- 
vided in Table 1. Note that each synaptic connection in the model 
is weighted by a connectivity constant denoted by Cspisp2 where 
SP1 and SP2, respectively, denote the source and target subpopu- 
lations. In addition, two Gaussian noise inputs pp{t)~ N{[ip, ap) 
and prc(f)~ N(u,tc, ore) were used to represent nonspecific 
inputs on pyramidal and thalamocortical cell subpopulations. 
Finally, signals simulated at the level of pyramidal cells in the 
cortical compartment and at the level of TC cells in the thalamic 
compartment were chosen as model outputs. They correspond 
to the sum of PSPs at each compartment respectively. The tem- 
poral dynamics of these signals provide a good approximation 
of actual LFPs. The model was implemented in Simulink®, 
and all other complementary scripts were implemented 
in MATLAB®. 

REAL DATA FOR MODEL TUNING 

We used real clinical data to tune the model into a func- 
tioning mode which simulates pathological activity. The clin- 
ical data set was limited to a unique patient who underwent 
thalamic DBS during the presurgical intracerebral EEG explo- 
ration (iEEG performed with depth electrodes implanted under 
stereotaxic conditions) at the Epilepsy Surgery Unit, Rennes 
University Hospital. This particular patient was chosen for 
two main reasons: (1) the pronounced frequency-dependent 
stimulation effects observed during his preoperative diagnos- 
tic iEEG exploration at LF, IF and HF in addition to (2) the 
existence of an epileptogenic zone in a limited area of the 
premotor cortex (PMC). 

In brief, this patient suffered from partial drug-resistant 
epilepsy since the age of two. MRI scans and EEG recordings 
pointed out the existence of a neuronal malformation known 
as FCD in the PMC at the origin of seizures. This type of 
cortical malformation is known for its epileptogenic features 
like neuronal hyperexcitability and hypersynchronization and its 
characteristic epileptiform discharges (continuous, rhythmic or 
semirhythmic spikes, and polyspikes) (Avoli et al., 2003; Palmini, 
2010) as shown in Figure 2C. Based on various clinical stud- 
ies reporting the modulation of epileptic cortical activity by the 
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Table 1 | Model parameters, values and interpretation. 



Parameter 


Value 


Interpretation 


Ac 


6 (optimized, pathological) 


Amplitude of the cortical average EPSP 




3 (normal) mV 




Ac 


K(t).A c mV 


Amplitude of the cortical average EPSP in response to thalamic input (only on 






subpopulation P) 


Be 


14 (optimized, pathological) 


Amplitude of the cortical average IPSP (GABA4 s j ow mediated currents) 




50 (normal) mV 




G c 


16.5 (optimized, 


Amplitude of the cortical average IPSP (GABA4 j ast mediated currents) 




pathological) 






zz tnormai/ mv 




A 

A Th 


0.0 mv 


Amplitude of the thalamic average EPSP 


n 

mh 


30 mV 


Amplitude of the thalamic average IPSP (GABA/\ s / ovv and GABAg receptors) 


Gjh 


22 mV 


Amplitude of the thalamic average IPSP (GABA4 j ast receptors) 


Am 


3.5 mV 


Amplitude of the reticular average EPSP 


x ac 


1/80s 


Time constant of cortical glutamate-mediated synaptic transmission. 


*bc 


1/35 s 


Time constant of cortical GABA-mediated synaptic transmission (GABA^ s i ow receptors) 


Igc 


1/180s 


Time constant of cortical GABA-mediated synaptic transmission (GABA^ f ast receptors) 


*aTh 


1/100s 


Time constant of thalamic glutamate-mediated synaptic transmission 


x bTh 


1/20s 


Time constant of thalamic GABA-mediated synaptic transmission (GABA^ s i ow and 






GABAg receptors) 


x gTh 


I / I DU S 


Time constant of thalamic GABA-mediated synaptic transmission (GABA4 f aS [ receptors) 




1/100 s 


Time constant of reticular glutamate-mediated synaptic transmission 




v 0 = 6mV, e 0 = 2.5 s _1 


Parameters of the nonlinear sigmoid function (transforming the average membrane 




r = 0.56mV-^ 


potential to an average density of action potentials) 


Cp-p 


135 


Collateral excitation connectivity constant 


Cp'-p 


108 


Collateral excitation connectivity constant 


c p-q 


33.75 


P to I2 connectivity constant 


c q-p 


33.75 


I2 to P connectivity constant 


c p-;f 


40.5 


P to connectivity constant 


C,c_,c 
'2 '1 


13.5 


/f to I2 connectivity constant 


c /f-p 


91.125 


to P connectivity constant 


Ctc-p 


120 


TC to P connectivity constant 


C TC-lC 


30 


TC to l9 connectivity constant 


C TC-f i 


45 


TC to I2 connectivity constant 


C TC-I? 


20 


TC to connectivity constant 


c rc-q< 


20 


TC to I2 connectivity constant 


C p _ l m 


30 


Pto /, connectivity constant 


C p ,Rt 
P-l 2 


30 


Pto connectivity constant 


Gp—TC 


20 


P to TC connectivity constant 


^l^-TC 


35 


1?* to TC connectivity constant 


Cm 


5 


1^ to TC connectivity constant 


M-P1 


U 


Mean of nonspecific cortical input 


\>-P2 


70 


Mean of nonspecific subcortical input 


apt 


20.v6 


Standard deviation of nonspecific cortical input 


Op? 


35.V6 


Standard deviation of nonspecific subcortical input 


Stc 


5 


Stimulation impact on subpopulation TC 


Smi 


4 


Stimulation impact on subpopulation /{" 


Spa 


4 


Stimulation impact on subpopulation 1™ 


fs 


1Hz- 150Hz 


Frequency of the stimulation signal (pulse train) 


Afs 


1 


Stimulation signal amplitude 



Model parameters used to reproduce LFPs FC o- Stimulation impact parameters S TC , S Btl and Sg t 2 are set to zero during the simulation of the NS scenario. These 
parameters are held constant for all other stimulation scenarios. 
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No 

Stimulation 
2 Hz 



3 
< 



Stimulation lAA^JWA^V^MnV 



50 Hz 
Stimulation 



70 Hz 
Stimulation 



100 Hz ^^|Vty\f''' ,lv 'w<V\^^ 
Stimulation 

Stimulation 



1 s 



FIGURE 2 | Frequency-dependent stimulation effects: real data. iEEG 
signals recorded during presurgical depth-EEG exploration in a patient with 
drug-resistant epilepsy. (A) MRI data showing the FCD (focal cortical 
dysplasia in the PMC) and the electrode trajectory. The red dot marks the 
position of the depth electrode in the FCD. (B) Zoom on the FCD. (C) DBS 
of the CMN modulated the pathological activity of the FCD in a 
frequency-dependent manner. LFS (2 Hz) and HFS (>70 Hz) suppressed 
pathological oscillations. IFS (50 Hz) had no effects. 



stimulation of the CM nucleus (Velasco et al., 1995, 1997, 2000, 
2001, 2007), itwas decided by neurologists and neursosurgeons to 
implant a depth electrode in this nucleus, as potentially beneficial 
for the patient who gave his informed consent. 

During the presurgical exploration, the stimulation of the tha- 
lamic CM nucleus (CMN) induced frequency-dependent modu- 
lation of the pathologic activity of the FCD (Figure 2). Readers 
may refer to (Pasnicu et al., 2013) for detailed information. 
Interestingly, LFS (2 Hz, 4 mA) and HFS (70, 100, and 150 Hz, 
0.8 mA) desynchronized the pathological activity of the FCD, 
while IFS (50 Hz, 0.8 mA) barely affected it. These segments of 



signals corresponding to either typical pathological activity or 
modulated activity (depending on stimulation conditions) were 
used to optimize the model parameters. 

PROCESSING OF REAL AND SIMULATED SIGNALS 

The use of signal processing techniques was necessary (i) to 
quantify the above-described effects of stimulation in real iEEG 
signals, and (ii) to define a feature-vector-based cost function for 
model parameter optimization. Figure 3A illustrates the feature 
extraction methodology. iEEG signals recorded in the FCD in 
absence of stimulation (LFPsfcd) and under different stimulation 
conditions were decomposed using an orthogonal matching pur- 
suit algorithm [matching pursuit toolkit — MPTK — (Krstulovic 
and Gribonval, 2006)]. First introduced in 1993 (Mallat and 
Zhifeng, 1993), matching pursuit is signal processing algorithm 
used to decompose any time series into a linear sum of wave- 
forms selected from a predefined dictionary based on a mother 
wavelet. To proceed, a proper multi-scalar dictionary of Gabor, 
Fourier, and Dirac atoms was first defined to account for real 
iEEG signal components (time-frequency atoms are waveforms 
well localized in both the time and the frequency domains). In line 
with (Krstulovic and Gribonval, 2006), the multi-scalar dictio- 
nary was formed by translation in time and amplitude/frequency 
modulation of atoms (defined as Gabor and Fourier functions 
in our case), over ten different user-defined time scales (i.e. the 
atom durations, ranging from 0.125 to 5 s). Then, the algorithm 
provided a table of time-frequency parameters associated to the 
detected atoms (atom type, central frequency, phase, scale, ampli- 
tude, position). Identified atoms were reconstructed using the 
extracted parameter table and their analytical expression. They 
were then associated to a given frequency band depending on 
their central frequency. These frequency bands corresponded to 
the classical EEG bands as defined in normal adults (8i [0 - 
1.9Hz],o 2 [1.9 -3.4 Hz], 0! [3.4 - 5.4 Hz], 6 2 [5.4 - 7.4Hz], «i 
[7.4 - 10Hz], a 2 [10 - 12Hz], Pi [12 -18 Hz], p 2 [18 - 24Hz], 
Y [24 - 128 Hz]) (Figure 3A, blue). Finally, a 9D feature vector 
Vp was defined from the normalized energy distribution in these 
frequency bands, itself computed as the sum of averaged (over 
time) atom energies relative to the total signal energy (Figure 3A, 
green). 

MODEL OPTIMIZATION UNDER THE "NO STIMULATION" CONDITION 

In order to simulate LFPsfcd, we optimized the excita- 
tion/inhibition ratio of the cortical compartment. Thus, the 
average EPSP/IPSP amplitude parameters of the cortical com- 
partment [Ac, Be, Gc} were considered as free parameters 
while all other model parameters were set to fixed values 
(Table 1). The optimization method is illustrated in Figure 3B. 
For each triplet {Ac, Be, Gc), the feature vector Vp, model of 
the model's output signal (cortical compartment's LFP) was cal- 
culated and compared to Vp rea i, i.e., the feature vector com- 
puted from the average of the 20 feature vectors, each com- 
puted on a 5s signal segment of real LFPsfcd- Feature vectors 
Vp, model and Vjj, re ai were computed as described in section 
Processing of Real and Simulated Signals. The optimization pro- 
cedure aimed at finding the triplet J Ac, Be, GcJ that min- 
imizes a cost function simply corresponding to the Euclidean 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



July 2013 | Volume 7 | Article 94 | 6 



Mina et al. 



Frequency-dependent effects of DBS 



iEEG signal 




MP 



Is 



Reconstructed signal 



Multi-scalar 
Dictionary 




8, 



Table of 
atom 
parameters 



Signal characterization 
and feature extraction 



Type 



Frequency f c 



Scale 



Amplitude 



Position 



Phase 



1s 



o 

> 

a) 



Normalized energy 
in frequency bands 



£(5 2 ) 
E(B,) 
£(9 2 ) 
BM 
E(a 2 ) 

£(Pi) 

e(W 

E(y) 



Association of 
atoms with the 
corresponding 
frequency bands 



ft, 



V F = 



V AA/ vAAAA/VW- 
'• wWVWVVAa- 



«2 



/S(£ ( ) 



Analytical 
reconstruction 
of atoms in the « 
corresponding 

frequency 
bands 



8, if f o e [0,1.9[ 
6, iff c e [1.9,3.4[ 
9~ if f c e[3.4,5.4[ 

iff c e[5.4,7.4[ 

o,iff c e[7.4,10[ 
a,iff c e[10,12.4[ 
Pi if f c e [12.4,181 
,if f c e [18,24[ 
y if f c G [24,128[ 



B 

Model 

parameter 

vector 




Ac 






B c 




Model 


[Gel 







Feature 
Extraction 



F, model 



t 



V, 



Freal 



} 



Euclidean 

distance 

minimization 



Parameter optimization 



FIGURE 3 | iEEG signal processing. (A) Feature vector extraction. Input 
signals were characterized using the matching pursuit (MP) method 
(dictionary of Gabor, Fourier, and Dirac atoms). Parameters of detected atoms 
(atom type, central frequency f c , scale, phase, amplitude, and position) are 
extracted by MP from input signals. Detected atoms are then associated with 
frequency bands (8! to y) depending on their proper central frequency. 
Sub-band (8! to y) signals were reconstructed from the sum of corresponding 



atoms, themselves obtained by fitting parameters into their analytic 
expression (see top left: input and reconstructed signals). The normalized 
energy vector [E(8i ) . . . £(y)]/(E(8i ) + ...+ E(y)] was chosen as the feature 
vector for further optimization of model parameters. (B) The model's free 
parameters Ac, Be, and Gc were optimized by minimizing the distance 
between the feature vector Vp mo dei of the simulated cortical LFP and the 
average of real feature vectors V F: regl of LFPs FCD . 



distance d(Vp, rea ;, Vp t mo dd) when parameters Ac, Be, and Gc 
span pre-defined ranges of values according to a Brute-Force 
procedure. 

RESULTS 

In this section, results regarding the identification of cellu- 
lar mechanisms underlying the modulation of cortical activity 
by thalamic DBS are reported. First, the model capability to 
reproduce signals similar to those recorded from the FCD in 
the patient was assessed, under two conditions (no stimula- 
tion and during stimulation). Three mechanisms contributing 
to frequency-dependant stimulation effects could be identified. 
Then, simulations were performed to analyze the marginal or 
joint contribution of these mechanisms at low, intermediate or 
high frequency stimulation. 



SIMULATION OF LFPs FCD UNDER NO STIMULATION CONDITION 

As a first step, we verified the ability of the model to gen- 
erate signals that resemble those recorded from the FCD in 
the considered patient {LFPs^qd)- This procedure, described in 
sections Processing of Real and Simulated Signals and Model 
Optimization Under the "No Stimulation" Condition, led us to 
identify a minimal distance (Figures 4A-C) and thus an optimal 

parameter vector |ac, Be, Gc J = {6, 14, 16.5} for which sim- 
ulated signals under the no stimulation condition have similar 
features as compared with those of real signals (Figure 4D). 

SIMULATION OF LFPs FCD UNDER STIMULATION CONDITIONS 

Actual LFPsfcd recorded at various stimulation frequencies (2, 
50, 70, 100, and 150 Hz) were first characterized using the 
matching pursuit method described in section Processing of 
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FIGURE 4 | Model parameter optimization. Normalized Euclidian 
distance between Vg real and Ve, model- Best fit (gray disk) between 
simulated and real LFPsfcd was obtained for (A) Ac = 6, (B) Be = 14, 
and (C) Gc = 16.5. (D) For these modified values of excitation and 



inhibition, the simulated signal exhibits similar characteristics as the iEEG 
signal recorded in the FCD. For standard values of excitation and 
inhibition {Ac = 3, Be = 50, Gc = 22), the model generates background 
EEG activity. 



Real and Simulated Signals (Figure 3A). Results are shown in 
Figure 5A where feature vectors of segments of LFPsfcd are 
represented in a 3D space where axes correspond to merged typ- 
ical EEG frequency bands (82 to 9i, 62 to Pi, P2 to y). Results 
show that the distribution of points in the 3D frequency space 
is not random but clustered, indicating that the frequency con- 
tent of LFPsfcd segments depends on the stimulation frequency. 
In addition, some clusters are very close. This is typically the 
case for i) the no stimulation (yellow) and the 50 Hz stimula- 
tion conditions (red) on the one hand, and ii) the 70 Hz (violet) 
and 150 Hz (cyan) stimulation conditions on the other hand. 
To go beyond the qualitative clustering performed by visual 
inspection of 3D plots, a K-means clustering algorithm imple- 
mented in MATLAB and using a Mahalanobis distance was used 
to automatically detect the three types of stimulation effects. 
Initial centroids were randomly chosen. The optimal cluster- 
ing that globally minimizes intra-cluster inertia is presented in 
Figure 5B. LFPsfcd segments were automatically classified into 
three subgroups. The first subgroup contains LFPsfcd segments 
obtained for low- frequency stimulation (LFS). The second sub- 
group gathers all segments recorded for high frequency stim- 
ulations (HFS, > 70 Hz). And finally, in the third subgroup, 
segments obtained under the no stimulation and the interme- 
diate stimulation frequency (IFS, 50 Hz) conditions are merged 



together, suggesting that this stimulation frequency does not 
reduce the "epileptiform aspect" of the activity reflected in the 
LFP. 

Based on this characterization of local field potentials recorded 
in the FCD (LFPsfcd), parameters Sjc-> Srh and Snt2 were man- 
ually tuned to lead the model to generate simulated signals which 
have spectral characteristics similar to those of actual LFPsfcd- 
Such a manual procedure was sufficient to reproduce stimulation 
effects observed in one patient. However, extending the study to 
a larger group of patients would have made imperative an auto- 
mated parameter fitting procedure based on the spectral charac- 
teristics of real EEG signals as in Rowe et al. (2004). Figure 5B 
shows the projection of representative simulated LFPsfcd in the 
3D frequency space ("M" triangles). As depicted, simulated sig- 
nals obtained for LFS, IFS and HFS were close to corresponding 
clusters obtained from real signals for the exact same computation 
of feature vectors. Shown in Figure 5C, these representatives sim- 
ulated LFPsfcd do not perfectly match actual signals. However, 
qualitatively similar bifurcations were observed in the model 
when the stimulation conditions are changed. Indeed, under the 
no stimulation (NS) and the IFS condition the model generates 
rhythmic slow oscillations (8) with superimposed faster activity 
(P, y), as observed in real data. For LFS and HFS conditions, 
strong modulation of this activity was also obtained in the model. 
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FIGURE 5 | Characterization and classification of real and simulated 
data. (A) 3-dimensional (3D) projection of feature vectors (V F rg3 j) 
corresponding to different stimulation conditions. This projection was 
obtained by summing some coordinates of initial 9D feature vectors to 
get 3D vectors [£(82)4-6(81), E(6 2 )+E(a 1 )+£(a 2 )+E(Pi ), 

F(fJ2)+£(y)l/[£(8i H h £(?)]■ Each vector was then represented by a 

point in the 3D space (82+81, 82 + 0! +02 + Pi, P2 + y)-Three main 
classes can be visually identified. (B) Clusters obtained using the 



k-means algorithm (Mahalanobis distance). The three clusters correspond 
to (i) low-frequency stimulation (LFS) effects (green squares), (ii) no 
stimulation (NS) and intermediate-frequency stimulation (IFS) effects 
(yellow squares), and (iii) high-frequency stimulation (HFS) effects (blue 
squares). Simulated signals corresponding to the four types of scenarios 
(NS, LFS, IFS, and HFS) were also projected in the same space 
(triangles). (C) Two-second segments of real and simulated signal during 
NS, LFS, IFS, and HFS. 



At LFS, in the model, the slow wave activity was strongly reduced 
but spike events occurred in the signals at the instant times of 
stimulation, mimicking, to some extent, comparable events also 
present in actual LFPsfcd- Finally, at HFS, slow oscillations (8) 
were abolished in the model which generates quasi-normal back- 
ground activity. This simulated activity was also comparable to 
real activity observed for HFS stimulation but disclosed less y 
activity. Note that these are the effects which were quantified in 
Figure 5B. The qualitative optimization procedure of parame- 
ters Stc> SRti, and Snt2 was then complemented by an evaluation 
of parameter sensitivity aimed at studying the impact of ran- 
dom changes affecting the parameter vector 0 = [Ac, Be, Gc, 
A-Th> Brh, Gjh, An t } on simulated signals. Parameter vector 0 
determines the excitability properties in the three model com- 
partments. As shown in Figure 6, results show that the simulated 
signals obtained under the four stimulation conditions (NS, LFS, 
IFS, HFS) stay "quite robust" (in the sense that waveforms are 

conserved) when parameters stay in the range [© ± t,. 0] with 
0 < t, < 0.2. 

MECHANISMS UNDERLYING FREQUENCY-DEPENDANT STIMULATION 
EFFECTS 

Three main mechanisms implemented in the model are required 
to mimic actually observed effects of the CM nucleus stimula- 
tion. These mechanisms are the following: (i) the presence of 
feed-forward inhibition (FFI) at the level of thalamic projections 



to the FCD, (ii) the presence of short-term depression (STD) at 
the level of the thalamocortical glutamatergic synapses and, (iii) 
the depolarization of RtN inhibitory interneurons targeting TC 
cells. 

This result raises an additional question: to what extent 
the joint effect of these mechanisms is necessary to reproduce 
frequency-dependant stimulation effects (LFS, IFS, and HFS). 
In order to assess their individual contribution, we performed 
simulations where each mechanism was either present in — or 
removed from — the model (the model parameters remaining 
unchanged). Results are displayed in Figure 7. First, they con- 
firmed that both FFI and STD mechanisms are jointly nec- 
essary in the model to suppress the epileptic activity in the 
FCD when LFS is being used since the withdrawal of either 
STD or FFI leads the model to generate epileptic activity at 
LFS. Second, results indicated that the RtN inhibitory interneu- 
rons targeting TC cells (both lf T and lf T subpopulations) 
must be affected (i.e., depolarized) by the stimulation to obtain 
a suppression of epileptic activity when HFS is being used, 
as observed in the patient. Third, and interestingly, an unex- 
pected effect was observed at IFS when the depolarization of 
if 1 interneurons was removed from the model. Indeed, epilep- 
tic activity was abolished in this case, which is really unlikely 
to occur during actual stimulation as both subtypes of neurons 
are expected to be affected by the direct stimulation of the CM 
nucleus. 
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These results were complemented by a deeper analysis of the 
thalamic output (i.e., the firing rate of TC cells) in response to 
stimulation at low, intermediate and high frequency. Results are 
provided in Figure 8. First, they showed that the thalamic out- 
put dramatically differs depending on the stimulation frequency 
(Figure 8A). Under the no stimulation condition, the firing rate 
continuously oscillates around a certain value (referred to as 
A, Figure 8A). At LFS, the firing rate was found to be lower, 
except at the stimulation times where it abruptly and transiently 
increased. At IFS, a balance was observed between time inter- 
vals for which the TC firing is above and below A. Finally, at 
HFS, the output of TC cells was found to be very low, i.e., sys- 
tematically under the threshold A. From these observations, we 
could define (i) two time intervals, Al and A2, for which the 
TC cells firing rate is either below A(A1) or above A (A2) 
and (ii) a "high to low firing" ratio (HtoLR) which provides 
an indication on the amount of time the TC cells spend fir- 
ing (up state) relatively to the amount of time they do not fire 
(down state). Figure 8B provides the evolution of the HtoLR 
when the stimulation frequency is progressively changing from 0 
to 150 Hz in the model. As depicted, these simulations indicated 
that three stimulation frequency ranges have dramatic effects on 
the firing of TC cells. First, from 0 to 20 Hz, the down state 



is predominant. Then, an abrupt jump was observed around 
22 Hz indicating that beyond this value, the firing rate dramat- 
ically increased. Interestingly, from 55 to 65 Hz, a progressive 
decrease of the HtoLR was observed. Then, after 70 Hz, the ratio 
is equal to zero indicating that TC cells did not fire anymore. 
Finally, in order to relate the thalamic activity with the corti- 
cal activity, we plotted the phase portraits (TC cell firing vs. 
cortical LFP) as illustrated in Figure 8C. Results confirmed the 
visual inspection of signals simulated at the two sites. For the 
no stimulation (NS) and for the intermediate frequency stimu- 
lation (IFS) conditions, phase portraits were found to be quite 
similar. They indicated the presence of mixed slow/fast oscilla- 
tions in both signals. For the low frequency stimulation (LFS) 
condition, oscillations in the simulated LFP in the FCD were 
reduced. They came along with short-duration, abrupt and rhyth- 
mic augmentations of the TC firing corresponding to stimulation 
pulses. Finally, for the high frequency stimulation (HFS) condi- 
tion, oscillations in both types of activity stayed confined to small 
amplitude values. 

DISCUSSION 

We modeled the thalamocortical loop in order to investigate 
frequency-dependent effects of electrical stimulation performed 
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FIGURE 7 I Conditions to reproduce frequency-dependent 
stimulation effects. Model output in the case where one of the 
implemented mechanisms (FFI, STD, depolarization of /™, and 
respectively) is removed at a time. LFS effects are not reproduced 
when the model does not account for FFI and STD. HFS effects 



in the thalamus and aimed at modulating the neocortical activ- 
ity. We chose to elaborate our model at a mesoscopic level, i.e., 
intermediate between microsocopic and macroscopic. 

Regarding the model architecture, we followed a similar 
approach to that used in previously proposed models of the tha- 
lamocortical loop (Robinson et al., 2002; Suffczynski et al., 2004; 
Breakspear et al., 2006; Roberts and Robinson, 2008; Marten et al., 
2009; Crunelli et al., 2011). Our model includes three main com- 
partments: cerebral cortex, reticular nucleus and thalamic relay. 
Subpopulations of neurons and interneurons located in these 
three structures interact via excitatory and/or inhibitory synaptic 
connections. The novelty with respect to aforementioned stud- 
ies is threefold. First, we modified the cortical compartment in 
order to better approximate the temporal dynamics of epileptic 
signals recorded in the FCD. This modification consisted in the 
use of two types of interneurons (mediating GABAergic IPSPs 
with slow and fast kinetics on cortical principal cells), as reported 
in a previous study (Molaee-Ardekani et al, 2010). Second, our 
model accounts for the direct effects of electrical stimulation. At 
this stage, we used the A V ~ \.E assumption according to which 
the perturbation of the mean membrane potential of neurons 
is a linear function of the electrical field magnitude induced by 
bipolar stimulation. This "X£" assumption was already used in 
neural mass models in the context of low-intensity direct hip- 
pocampal stimulation to anticipate seizures (Suffczynski et al., 
2008) as well as in the analysis of the stimulus-response relation- 
ship of DBS in healthy animals (Adhikari et al., 2009). However, 
it is worth mentioning that in our model, the three subtypes of 
neurons (TC cells and both subpopulations of inhibitory neurons 
in the RtN) are depolarized by the stimulation, as suggested in 



require the depolarization of both reticular populations /j and /"'. 
Suppression of epileptic activity is observed at IFS when ff' 
interneurons are removed. Red dotted lines indicate situations where 
simulated signals do not match real ones for given stimulation 
condition. 



(Molaee-Ardekani et al., 2013) and conversely to (Adhikari et al., 
2009) where only principal cells are impacted. And third, our 
model includes two well-known mechanisms at the cortical level: 
feed-forward inhibition (FFI) and short-term depression (STD). 

As in any modeling approach, our approach has some limita- 
tions. First, the chosen modeling level does not allow for analyzing 
sub-cellular mechanisms involved in stimulation-evoked changes. 
Similarly, it does not account for direct activation of axons by 
stimulation versus somatic inhibition (Mclntyre et al., 2004b), 
nor for the mechanisms of orthodromic/antidromic propaga- 
tion of action potentials due to stimulation (Degos et al., 2005; 
Hammond et al., 2007; Dorval et al., 2008). Second, a strong 
assumption in the type of model we used (neural mass) is related 
to the intrinsic synchronization among neurons included in a 
given sub-population. This assumption does not allow for rep- 
resenting either de- or weakly-synchronized firing patterns that 
may be observed during epileptic activity, in particular during 
high frequency oscillations that can be encountered in FCDs 
(Brazdil et al., 2010). Nevertheless, we could accurately repro- 
duce the abnormal rhythms generated in the FCD suggesting that 
main pyramidal cells have a relatively synchronized activity in 
this epileptogenic tissue. Third, regarding plasticity- related mech- 
anisms, we only implemented short-term effects (i.e., STD) and 
neglected long-term plastic changes that may be induced by DBS 
(Shukla etal, 2013). 

Despite these limitations, we could identify a number of 
mesoscopic factors which could explain the frequency-dependent 
mechanisms of thalamic stimulation. The model was tuned using 
electrophysiological data recorded in a patient in whom the cen- 
tromedian nucleus (CMN) stimulation was particularly efficient 
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to reduce the epileptic activity of a FCD located in the premo- 
tor cortex, in a frequency-specific manner. The main findings are 
summarized in Figure 9. 

"NO STIMULATION" (NS) CONDITION 

In the model, under the NS condition, excitation among pyra- 
midal cells had to be increased and inhibition had to be reduced 
in the cortical compartment for producing "pathological" oscilla- 
tory rhythms, as observed in the FCD. The thalamocortical loop 
was found to be responsible for these pathological dynamics, 
characteristic of FCDs. These findings are in line with histologi- 
cal studies showing that these typical oscillations are generated in 
altered brain tissue, where inhibition is partially deteriorated or 
dysfunctioning (Calcagnotto et al., 2005), and where excitation is 
heavily increased (Avoli et al., 2003). In addition to neuron alter- 
ations in the dysplastic tissue (Sisodiya et al., 2009), FCD keeps 
sufficient projections to — and input from — other brain struc- 
tures to propagate pathological dynamics (Avoli et al., 2003). As 
mentioned, the presence of connections with subcortical struc- 
tures was a necessary condition in the model for producing 



pathological oscillations resembling those actually recorded in the 
FCD (Figure 9A). 

LOW-FREQUENCY STIMULATION (LFS) CONDITION 

For the low-frequency stimulation (LFS, / < 20 Hz) condition, 
two mechanisms were found to play a major role for the abor- 
tion of epileptic activity in the FCD: short-term depression (STD, 
i.e., decreased excitatory synaptic efficacy in thalamus-to-cortex 
connections) and feed-forward inhibition (FFI, i.e., excitation of 
inhibitory cortical interneurons by TC cells) (Figure 9B). 

STD was reported in previous studies concerning cortical 
adaptation to thalamic stimulation, and suggesting that electrical 
LFS of TC cell axons in vivo resulted in a 40% reduction in cor- 
tical EPSPs (Chung et al., 2002). In the same context, LFS trains 
in adult anaesthetized rats provoked transient long-term depres- 
sion of thalamocortical synapses; this was measured by up to 40% 
drop in cortical EPSPs after LFS trains and under the effect of 
GABA antagonist (Speechley et al, 2007). 

As mentioned above, the LFS effects could not be repro- 
duced by the model without incorporating also FFI. Actually, 
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FIGURE 9 | Frequency-dependent mechanisms underlying DBS. (A) 

Under the no stimulation (NS) condition, the thalamocortical loop is 
responsible for pathological oscillatory rhythms observed in the FCD. (B) For 
low-frequency stimulation (LFS), feed-forward inhibition (FFI, i.e., excitation of 
inhibitory cortical interneurons by TC cells) and short-term depression (STD, 
i.e., decreased excitatory synaptic efficacy in thalamus-to-cortex connections) 
was found to play a major role for the abortion of epileptic activity in the FCD. 



(C) For the intermediate-frequency stimulation (IFS) condition, thalamic 
output is reinforced (increase of TC cells firing) leading to an increase of the 
average excitatory post-synaptic potential (EPSP) on cortical pyramidal cells 
and to no "anti-epileptic" effect. (D) For high-frequency stimulation (HFS), the 
direct and sustained excitation of reticular nucleus (RtN) interneurons leads 
to dramatic decrease of TC cells firing rate and to a suppression of epileptic 
activity. 



thalamocortical ascending fibers directly target pyramidal neu- 
rons as well as cortical GABAergic interneurons inducing EPSPs 
in both cell types (Pouille and Scanziani, 2001). In the model, 
while less efficient (STD) thalamic EPSPs arrive directly onto 
pyramidal neurons, IPSPs induced by thalamic stimulation also 
arrive on pyramidal neurons (FFI) lagging by 1-2 ms. This short 
latency between the onset of thalamocortical excitation and 
the onset of feed-forward inhibition presents a temporal "win- 
dow of opportunity" for pyramidal cells to integrate excitatory 
and inhibitory inputs, thus keeping the transmembrane poten- 
tial below firing threshold. In the literature, neuroanatomical 
and neurophysiological studies (Isaacson and Scanziani, 2011) 
showed the functional importance of FFI in regulating corti- 
cal dynamics by controlling cortical excitability (Cabernet et al., 
2005). Our study suggests that LFS regulates cortical excitability 
by a dual mechanism of FFI and STD (Figure 9B). 

INTERMEDIATE-FREQUENCY STIMULATION (IFS) CONDITION 

For the intermediate-frequency stimulation (IFS, 20 < 
/ < 70 Hz) condition, results indicated that the thalamic 
output is reinforced (increase of TC cells firing) and leads to an 
increase of the average excitatory post-synaptic potential (EPSP) 
on cortical pyramidal cells (Figure 9C). This effect corresponds 
to an increase of the spatiotemporal summation of unitary 
EPSPs. In this case, both the cortical excitability and the gain in 
the excitatory thalamocortical loop is increased, leading to "no 
anti-epileptic" effect. We did not find much studies using DBS 
stimulation in the intermediate frequency range of (20-60 Hz) in 
the context of epilepsy. Nevertheless, it is noteworthy that 50 Hz 
stimulation frequency is classically used during the presurgical 
evaluation of patient with intractable partial epilepsy in order to 



trigger seizures and delineate the epileptogenic zone (Talairach 
et al., 1974; Jayakar et al, 1992). The same frequency range is also 
known to provoke afterdischarges and was actually used in the 
kindling model of epilepsy (Goddard, 1967; Racine, 1972). 

HIGH-FREQUENCY STIMULATION (HFS) CONDITION 

Finally, for the high-frequency stimulation (HFS, / > 70 Hz) con- 
dition, the direct and sustained excitation of reticular nucleus 
(RtN) interneurons leads to strong inhibition of TC cells and 
thus to dramatic decrease of their firing rate. Despite the fact 
that TC neurons are also affected by stimulation, the response 
of reticular GABAergic neurons to stimulation and the higher 
efficiency of GABA-mediated currents ensure that IPSPs over- 
ride EPSPs on TC cells. In this case, the reduced excitatory 
input to cortical pyramidal cells also leads to a suppression of 
epileptic activity (Figure 9D). This result corroborates reported 
stimulation studies where HFS (>100Hz) was associated with 
significant decrease in epileptiform discharges in vitro, and reduc- 
tion in seizure frequency in responding patients (Velasco et al., 
2006; Fisher et al., 2010). This hypothesis is in line with recent 
findings suggesting that HFS of the globus pallidus (GPi) in dys- 
tonia patients decreased its firing by stimulation-evoked GABA 
release from afferent fibers and thereby the enhancement of 
inhibitory synaptic transmission by HFS (Liu et al., 2012). 
Similarly, HFS (100Hz-130Hz) of the STN neurons in vitro 
showed a suppression of the activity of the majority of neu- 
rons by the reinforcement of inhibitory responses (Filali et al., 
2004). Other HFS studies also provided evidence on the inhibi- 
tion of GPi output during HFS in human patients (Dostrovsky 
et al., 2000) as well as the disruption thalamocortical network's 
dysrhythmia (Mclntyre and Hahn, 2010; Kendall et al, 2011). 
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CONCLUSION 

In epilepsy research, it is well-admitted that there is, unfor- 
tunately, a lack of tangible results regarding the effects of 
electrical stimulation in the brain. Therefore, the very crucial 
issue of choosing the "optimal" stimulation parameters remains 
unsolved, whatever the stimulation procedure. Although com- 
putational models are always based on a number of simplifying 
assumptions, we think that they provide an efficient framework to 
(i) account for the many and essential factors that may intervene 
during stimulation procedures and (ii) analyze the links between 
these factors in a formal manner. This approach is particularly 
fruitful when models are well grounded in experimental/clinical 
data (Wendling et al, 2012). This is somehow a weak point of this 
study since we could make use of data sets recorded in one patient 
only. However, it should be mentioned that these very informa- 
tive data sets stay relatively rare since many conditions have to be 
met (patient candidate to surgery, FCD, electrodes positioned in 
appropriate structures). 

At this stage, the face value of the model is satisfactory. The 
next step is obviously to test the model predictions using animal 



models. Experiments can be undertaken in rodents with elec- 
trodes implanted in the cerebral cortex and in the thalamus. 
First, we could start with control animals to assess the modula- 
tion of cortical rhythms during/after direct thalamic stimulation 
at various frequencies and for controlled vigilance states (sleep, 
awake, resting, exploratory). In these controls, some drugs can 
be used to alter some parameters related to synaptic transmission 
(in a more or less specific manner) which have a correspondence 
in the model, on the other hand. Then, refined experimental 
models could be introduced to get closer to the epilepsy context 
including models of developmental dysplastic lesions [see review 
in Schwartzkroin and Wenzel (2012)]. Hopefully, this combined 
computational/experimental approach will help us to disclose 
some of the highly intricate effects of DBS either at local or at 
network level. 
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